#########################################################
## Install and load missing R packages [ packages(x) ]
## Written by TszKin Julian
## at https://rlearner.wordpress.com
#########################################################

packages<-function(x){
	   x<-as.character(match.call()[[2]])
	   if (!require(x,character.only=TRUE)){
      install.packages(pkgs=x,repos="http://cran.r-project.org")
      require(x,character.only=TRUE)}}


#########################################################
## Convert meters to degrees [ m2d(m) ]
#########################################################


m2d <- function(m){
out <- (m/1852)/60
return(out)
}
deg <- m2d(500000) ## 500 km

#########################################################
## Row normalization [ row.norm(W) ]
## * The following function accepts a binary matrix (W)
## * It produces a row-standardized weight matrix
#########################################################

row.norm <- function(W){
	for(k in 1:nrow(W)){
	if(sum(W[k,]>0))
	W[k,] <- W[k,]/sum(W[k,])
	}
	return(W)
	}

#########################################################
## Standard deviation above mean [ sdmean(x) ]
## * The following function accepts a numeric vector (x)
## * It produces a value 1 Std.Dev. above the mean
#########################################################

sdmean <- function(x) mean(x,na.rm=T) + sd(x,na.rm=T)

#########################################################
## Standard deviation shift [ sdshift(modelobj,shift) ]
## * The following function accepts a Zelig model
## * object based on Model 1/Table 1 from  
## * GW2006 (modelobj), a desired value shift (shift)
## * It produces a matrix of changes in predicted 
## * transition probabilities resulting from the shift
## * and 95% CI's
#########################################################

sdshift <- function(mod, shift) {

## Extract names of key terms

v1 <- attr(mod$terms,"term.labels")[2]
v2 <- attr(mod$terms,"term.labels")[9]
v2. <- v1
v3 <- attr(mod$terms,"term.labels")[14]
v3. <- gsub("I(","",v3,fixed=T)
v3. <- gsub(" * bautlag)","",v3.,fixed=T)


## Calculate SD shifts

shift <- apply(mod$data,2,sdmean)

## Set baseline values for simulation
## * assumptions: no civil war, no neighboring transitions

baseBeta <- setx(mod, bautlag=0, cwar=0) 
baseGamma <- setx(mod, bautlag=1, cwar=0)
baseBeta[v3] <- 0
baseGamma[v3] <- 0

s.Beta <- sim(mod, baseBeta, num=10000)
s.Gamma <- sim(mod, baseGamma, num=10000)

## Counterfactual 1: proportion of neighboring democracies

cf1Beta <- setx(mod, bautlag=0, cwar=0) 
cf1Gamma <- setx(mod, bautlag=1, cwar=0)
cf1Beta[v3] <- 0
cf1Gamma[v3] <- 0

cf1Beta[v1] <- shift[v1]
cf1Gamma[v1] <- shift[v1]

cf1Beta[v2] <- shift[v2.]*0
cf1Gamma[v2] <- shift[v2.]*1

cf1Beta. <- sim(mod, cf1Beta, num=10000)
cf1DA <- c(mean(cf1Beta.$qi$ev-s.Beta$qi$ev), quantile(cf1Beta.$qi$ev-s.Beta$qi$ev, .025), quantile(cf1Beta.$qi$ev-s.Beta$qi$ev, .975))

cf1Gamma. <- sim(mod, cf1Gamma, num=10000)
cf1AD <- c(mean((1-cf1Gamma.$qi$ev)-(1-s.Gamma$qi$ev)), quantile((1-cf1Gamma.$qi$ev)-(1-s.Gamma$qi$ev),.025), quantile((1-cf1Gamma.$qi$ev)-(1-s.Gamma$qi$ev),.975))

## Counterfactual 2: neighboring transition

cf2Beta <- setx(mod, bautlag=0, cwar=0) 
cf2Gamma <- setx(mod, bautlag=1, cwar=0)

cf2Beta[v3] <- shift[v3.]*0
cf2Gamma[v3] <- shift[v3.]*1

cf2Beta. <- sim(mod, cf2Beta, num=10000)
cf2DA <- c(mean(cf2Beta.$qi$ev-s.Beta$qi$ev), quantile(cf2Beta.$qi$ev-s.Beta$qi$ev, .025), quantile(cf2Beta.$qi$ev-s.Beta$qi$ev, .975))

cf2Gamma. <- sim(mod, cf2Gamma, num=10000)
cf2AD <- c(mean((1-cf2Gamma.$qi$ev)-(1-s.Gamma$qi$ev)), quantile((1-cf2Gamma.$qi$ev)-(1-s.Gamma$qi$ev),.025), quantile((1-cf2Gamma.$qi$ev)-(1-s.Gamma$qi$ev),.975))

## Output

results <- rbind(cf1DA, cf1AD, cf2DA, cf2AD)
colnames(results) <- c("Mean","Lower","Upper")
return(results)

}

#########################################################
## Area under the ROC curve [ AUC(y,fitted) ]
## * The following function calculates the area under
## * the Receiver-Operator Characteristic (ROC) curve.
## * It accepts as an argument a vector of observed 
## * values (y) and fitted values (fitted), which must be
## * of same length and order.
## * It generates an AUC statistic, which ranges from
## * 0 (worst) to .5 (random) to 1 (perfect fit).
#########################################################


AUC <- function(y, fitted){
    n1 <- sum(y)
    n <- length(y)
    out <- (mean(rank(fitted)[y == 1]) - (n1 + 1)/2)/(n - n1)
    return(out)
}

##########################################################
## Model Fitting Function [ models(X) ]
## * The following function fits the GW2006 
## * model for each of the 14 neighbor definitions
## * It accepts as an argument an in-sample dataset (X).
## * It generates a list of 14 Zelig model objects, which 
## * can be used for prediction and simulation
##########################################################

 

models <- function(X){
	
## Geographic: contiguity

m0.cont <- zelig(I(baut) ~ log(laggdppc) + pnbdem 
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(pnbdem*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(nbtd*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)

## Geographic: minimum distance
 
m0.mdn <- zelig(I(baut) ~ log(laggdppc) + Wpropdem.dist.2 
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(Wpropdem.dist.2*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(Wdemtr.dist.2*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)

## Geographic: K=4 nearest neighbors
 
m0.k4 <- zelig(I(baut) ~ log(laggdppc) + Wpropdem.k4.2 
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(Wpropdem.k4.2*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(Wdemtr.k4.2*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)

## Geographic: sphere of influence

m0.soi <- zelig(I(baut) ~ log(laggdppc) + Wpropdem.soi.2 
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(Wpropdem.soi.2*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(Wdemtr.soi.2*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)

## Network: trade (minumum distance)
 
m1.tr.mdn <- zelig(I(baut) ~ log(laggdppc) + W.pdem.mdn.tr
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.mdn.tr*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.mdn.tr*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)

## Network: trade (K=4 nearest neighbors)
 
m1.tr.k4 <- zelig(I(baut) ~ log(laggdppc) + W.pdem.k4.tr
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.k4.tr*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.k4.tr*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)

## Network: trade (pseudo sphere of influence)
 
m1.tr.soi <- zelig(I(baut) ~ log(laggdppc) + W.pdem.soi.tr
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.soi.tr*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.soi.tr*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)
 
## Network: ethnicity (minumum distance)

m1.eth.mdn <- zelig(I(baut) ~ log(laggdppc) + W.pdem.mdn.eth
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.mdn.eth*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.mdn.eth*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)
 
## Network: ethnicity (K=4 nearest neighbors)
 
m1.eth.k4 <- zelig(I(baut) ~ log(laggdppc) + W.pdem.k4.eth
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.k4.eth*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.k4.eth*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)
 
## Network: ethnicity (pseudo sphere of influence)
 
m1.eth.soi <- zelig(I(baut) ~ log(laggdppc) + W.pdem.soi.eth
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.soi.eth*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.soi.eth*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)
 
## Network: IGO (minumum distance)

m1.igo.mdn <- zelig(I(baut) ~ log(laggdppc) + W.pdem.mdn.igo
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.mdn.igo*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.mdn.igo*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)
 
## Network: IGO (K=4 nearest neighbors)
 
m1.igo.k4 <- zelig(I(baut) ~ log(laggdppc) + W.pdem.k4.igo
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.k4.igo*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.k4.igo*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)
 
## Network: IGO (pseudo sphere of influence)
 
m1.igo.soi <- zelig(I(baut) ~ log(laggdppc) + W.pdem.soi.igo
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.soi.igo*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.soi.igo*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)
 
## Network: Alliance (binary)
 
m1.al <- zelig(I(baut) ~ log(laggdppc) + W.pdem.al
  + cwar + ipyears + egr + propdem 
  + I(bautlag) + I(log(laggdppc)*bautlag) 
  + I(W.pdem.al*bautlag) + I(cwar*bautlag) + I(ipyears*bautlag) 
  + I(egr*bautlag) + I(propdem*bautlag)
  + I(W.demtr.al*bautlag), model="probit", data=X,
 trace=1, maxit=10000, cite=F)
 
return(list(m0.cont,m0.mdn,m0.k4,m0.soi,m1.tr.mdn,m1.tr.k4,m1.tr.soi,m1.eth.mdn,m1.eth.k4,m1.eth.soi,m1.igo.mdn,m1.igo.k4,m1.igo.soi,m1.al))
}



##########################################################
## Regime Change Prediction [ regimepred(m0.data,modz,mats) ]
## * The following function accepts scenario data (m0.data),
## * a list of model objects (modz) and
## * a list of weight matrices (mats).
## * It generates a matrix of predicted probabilities
##########################################################
i <- 12
regimepred <- function(modz,m0.data,mats,labels){

preds.mean <- matrix(NA,nrow=nrow(m0.data),ncol=15)
preds.mean[,1] <- m0.data$CNTRY_CODE
colnames(preds.mean) <- c("CNTRY_CODE",labels)
preds.l <- preds.mean
preds.u <- preds.mean

for(i in 1:14){
print(i)
m0 <- modz[[i]]
# Extract relevant variable names
v1 <- attr(modz[[i]]$terms,"term.labels")[2]
v2 <- attr(modz[[i]]$terms,"term.labels")[9]
v2. <- v1
v3 <- attr(modz[[i]]$terms,"term.labels")[14]
v3. <- gsub("I(","",v3,fixed=T)
v3. <- gsub(" * bautlag)","",v3.,fixed=T)

if(mean(m0.data$numid==rownames(mats[[i]]))==1){

	
# Prop of neighboring democracies
m0.data[,v1] <- 1-row.norm(mats[[i]])%*%m0.data$baut

# Number of neighboring transitions
m0.data[,v3.] <- mats[[i]]%*%(m0.data$baut==0 & m0.data$bautlag==1)
}

# Create model matrix

X <- cbind(1,log(m0.data$laggdppc),m0.data[,c(v1,"cwar",
"ipyears","egr","propdem","bautlag")],
m0.data$bautlag*log(m0.data$laggdppc),
m0.data$bautlag*m0.data[,v1],
m0.data$bautlag*m0.data$cwar,m0.data$bautlag*m0.data$ipyears,
m0.data$bautlag*m0.data$egr,m0.data$bautlag*m0.data$propdem,
m0.data$bautlag*m0.data[,v3.])
colnames(X) <- c("(Intercept)",attr(m0$terms,"term.labels"))
X <- as.matrix(X)
for(k in 1:15){
	X[,k] <- as.numeric(X[,k])
	}
	
# Run simulation

Y.t1 <- m0.data$bautlag
M <- modz[[i]]
B <- M$coefficients
VCOV <- summary(M)$cov.scaled
pred <- matrix(NA, nrow=nrow(X),ncol=1000)
set.seed(12345)

for(m in 1:1000){
beta <- rmvnorm(n=1,mean=B[1:7],sigma=VCOV[1:7,1:7])
alpha <- rmvnorm(n=1,mean=B[8:14],sigma=VCOV[8:14,8:14])
nbtd <- rnorm(n=1,mean=B[15]/2,sd=(VCOV[15,15])^2/2)
coef.AA <- c(beta+alpha,nbtd)
coef.DA <- c(beta,nbtd)
pred[,m] <- ifelse(Y.t1==1,dnorm(X[,c(1:7,15)]%*%coef.AA),dnorm(X[,c(1:7,15)]%*%coef.DA))
#pred[,m] <- dnorm(X%*%t(rmvnorm(n=1,mean=B,sigma=VCOV)))
}

# Compute summary statistics

preds.mean[,i+1] <- apply(X=pred,MARGIN=1,FUN=mean)
preds.l[,i+1] <- apply(X=pred,MARGIN=1,FUN=function(x){quantile(x,.025)})
preds.u[,i+1] <- apply(X=pred,MARGIN=1,FUN=function(x){quantile(x,.975)})
}

# Return matrices

return(list(preds.mean,preds.l,preds.u))
}




##########################################################
## Regime Change Simulation [ regimesim(cnt,modz,mats,labels) ]
## * The following function accepts country names (cnt)
## * for the locations of regime change,
## * a list of model objects (modz) and
## * a list of weight matrices (mats).
## * It generates a list of first differences, 
## * predicted probabilities and
## * confidence intervals for each type of transition.
##########################################################



regimesim <- function(cnt,modz,mats,labels=c(1:length(mats))){



# Create empty list of transition probabilities

preds <- list()


## Extract a single cross-section of the data, clean it up

m0 <- modz[[1]]
m0.data <- subset(m0$data, m0$data$year == 1998)
m0.data <- subset(m0.data, m0.data$baut != "NA" | m0.data$bautlag != "NA" | m0.data$demtr != "NA")
m0.data <- m0.data[m0.data[,"numid"]%in%colnames(mats[[1]]),]

# Find row corresponding to country(ies)

cntnum <- which(m0.data$MAP.CNTRY%in%cnt)


# Set regime type (AA)

m0.data[cntnum,"baut"] <- 1
m0.data[cntnum,"bautlag"] <- 1
m0.data[-cntnum,"baut"] <- m0.data[-cntnum,"bautlag"]
m0.data$demtr <- ifelse(m0.data$baut==0 & m0.data$bautlag==1,1,0)


# Compute Probabilities

pred <- regimepred(modz,m0.data,mats,labels)
preds[["AA"]] <- pred[[1]]
preds[["AA.l"]] <- pred[[2]]
preds[["AA.u"]] <- pred[[3]]


# Set regime type (AD)

## Extract a single cross-section of the data, clean it up

m0 <- modz[[1]]
m0.data <- subset(m0$data, m0$data$year == 1998)
m0.data <- subset(m0.data, m0.data$baut != "NA" | m0.data$bautlag != "NA" | m0.data$demtr != "NA")
m0.data <- m0.data[m0.data[,"numid"]%in%colnames(mats[[1]]),]

m0.data[cntnum,"baut"] <- 0
m0.data[cntnum,"bautlag"] <- 1
m0.data[-cntnum,"baut"] <- m0.data[-cntnum,"bautlag"]
m0.data$demtr <- ifelse(m0.data$baut==0 & m0.data$bautlag==1,1,0)

# Compute Probabilities

pred <- regimepred(modz,m0.data,mats,labels)
preds[["AD"]] <- pred[[1]]
preds[["AD.l"]] <- pred[[2]]
preds[["AD.u"]] <- pred[[3]]


# Set regime type (DA)

## Extract a single cross-section of the data, clean it up

m0 <- modz[[1]]
m0.data <- subset(m0$data, m0$data$year == 1998)
m0.data <- subset(m0.data, m0.data$baut != "NA" | m0.data$bautlag != "NA" | m0.data$demtr != "NA")
m0.data <- m0.data[m0.data[,"numid"]%in%colnames(mats[[1]]),]

m0.data[cntnum,"baut"] <- 1
m0.data[cntnum,"bautlag"] <- 0
m0.data[-cntnum,"baut"] <- m0.data[-cntnum,"bautlag"]
m0.data$demtr <- ifelse(m0.data$baut==0 & m0.data$bautlag==1,1,0)

# Compute Probabilities

pred <- regimepred(modz,m0.data,mats,labels)
preds[["DA"]] <- pred[[1]]
preds[["DA.l"]] <- pred[[2]]
preds[["DA.u"]] <- pred[[3]]


# Set regime type (DD)

## Extract a single cross-section of the data, clean it up

m0 <- modz[[1]]
m0.data <- subset(m0$data, m0$data$year == 1998)
m0.data <- subset(m0.data, m0.data$baut != "NA" | m0.data$bautlag != "NA" | m0.data$demtr != "NA")
m0.data <- m0.data[m0.data[,"numid"]%in%colnames(mats[[1]]),]

m0.data[cntnum,"baut"] <- 0
m0.data[cntnum,"bautlag"] <- 0
m0.data[-cntnum,"baut"] <- m0.data[-cntnum,"bautlag"]
m0.data$demtr <- ifelse(m0.data$baut==0 & m0.data$bautlag==1,1,0)

# Compute Probabilities

pred <- regimepred(modz,m0.data,mats,labels)
preds[["DD"]] <- pred[[1]]
preds[["DD.l"]] <- pred[[2]]
preds[["DD.u"]] <- pred[[3]]


# Calculate tables of first differences


m0 <- modz[[1]]
m0.data <- subset(m0$data, m0$data$year == 1998)
m0.data <- subset(m0.data, m0.data$baut != "NA" | m0.data$bautlag != "NA" | m0.data$demtr != "NA")
m0.data <- m0.data[m0.data[,"numid"]%in%colnames(mats[[1]]),]


vizs <- lapply(2:15,function(wnum){
viz <- matrix(NA, nrow=nrow(m0.data), ncol=13)
colnames(viz) <- c("CNTRY_CODE","CNTRY_NAME","R_AD","R_ADl","R_ADu","R_DA","R_DAl","R_DAu","bautlag","preds.AA","preds.AD","preds.DD","preds.DA")
viz[,1] <- m0.data$CNTRY_CODE
viz[,2] <- as.character(m0.data$CNTRY_NAME)
viz[,3] <- preds[["AD"]][,wnum] - preds[["AA"]][,wnum] ## Country Democratizes
viz[,4] <- preds[["AD.l"]][,wnum] - preds[["AA.l"]][,wnum] 
viz[,5] <- preds[["AD.u"]][,wnum] - preds[["AA.u"]][,wnum] 
ci.l <- apply(matrix(as.numeric(as.character(viz[,4:5])),ncol=2),1,min)
ci.u <- apply(matrix(as.numeric(as.character(viz[,4:5])),ncol=2),1,max)
viz[,4] <- ci.l
viz[,5] <- ci.u
viz[,6] <- preds[["DA"]][,wnum] - preds[["DD"]][,wnum] ## Country Autocratizes
viz[,7] <- preds[["DA.l"]][,wnum] - preds[["DD.l"]][,wnum]
viz[,8] <- preds[["DA.u"]][,wnum] - preds[["DD.u"]][,wnum]
ci.l <- apply(matrix(as.numeric(as.character(viz[,7:8])),ncol=2),1,min)
ci.u <- apply(matrix(as.numeric(as.character(viz[,7:8])),ncol=2),1,max)
viz[,7] <- ci.l
viz[,8] <- ci.u
viz[,9] <- m0.data$bautlag
viz[,10] <- preds[["AA"]][,wnum]
viz[,11] <- preds[["AD"]][,wnum]
viz[,12] <- preds[["DD"]][,wnum]
viz[,13] <- preds[["DA"]][,wnum]

viz <- as.data.frame(viz)
for(i in 3:9){
	viz[,i] <- as.numeric(as.character(viz[,i]))
	}
viz})
names(vizs) <- labels


# Return list

return(vizs)
}


##########################################################
## Likelihood function for SAR model [ sar.lik(theta,y,X,W) ] 
##########################################################


  sar.lik <-function(theta,y,X,W){
  	sigma2 <- theta[1]^2
		rho <- theta[2]
		betas <- theta[3:length(theta)]
		n <- nrow(X)
		eps <- (diag(n) - rho*W)%*%y - X%*%betas
		J <- sum(log(1-rho*eigen(W)$values))
		lnL <- -(n/2)*log(pi*sigma2)+J-(t(eps)%*%eps)/(2*sigma2)
		return(-lnL)
	}

##########################################################
## MLE for SAR [ my.sar(y,X,W) ]
## * Maximum likelihood estimation for a SAR model 
## * Takes as an input a vector of outcomes (y),
## * a matrix of predictor variables (X),
## * and a spatial weights matrix (W).
## * Output is a list of parameters and summary statistics.
##########################################################

my.sar <- function(y,X,W){
	init <- c(.1,.1,rep(.1,ncol(X)))
	lb <- c(1e-20,-.999,rep(-1e10,ncol(X)))
	ub <- c(1e20,.999,rep(1e10,ncol(X)))
#pars <- optim(init,y=y,X=X,W=W1,fn=sar.lik,method="L-BFGS-B",lower=lb,upper=ub,hessian=T,control=list(parscale=c(.1,.1,rep(.1,ncol(X)))))
pars <- nlminb(start=init, objective=sar.lik, gradient = NULL, hessian = sar.lik, y,X,W, lower = lb, upper = ub)
hes <- hessian(func=sar.lik,x=pars$par,y=y,X=X,W=W)
s2.hat <- pars$par[1]
rho.hat <- pars$par[2]
beta.hat <- pars$par[3:length(pars$par)]
LL <- -pars$objective
aic <- 2*ncol(X)-2*LL
resid <- y-(rho.hat*W%*%y+X%*%beta.hat)
SSR <- sum(resid^2)
OI <- solve(hes)
se <- sqrt(diag(OI))
coef <- cbind(c(rho.hat,beta.hat),se[2:length(se)])
#coef <- round(coef,12)
colnames(coef) <- c("Coefficient","Standard Error")
rownames(coef) <- c("Rho",paste("X",1:ncol(X),sep=""))
return(list(pars=pars,coef=coef,s2=s2.hat,LL=LL,AIC=aic,SSR=SSR))
}

##########################################################
## Likelihood function for SAR w/ 2 W's [ sar2.lik(theta,y,X,W1,W2) ] 
##########################################################


  sar2.lik <-function(theta,y,X,W1,W2){
    sigma2 <- theta[1]^2
		rho1 <- theta[2]
		rho2 <- theta[3]
		betas <- theta[4:length(theta)]
		n <- nrow(X)
		eps <- (diag(n) - rho1*W1 - rho2*W2)%*%y - X%*%betas
		J <- sum(log(1-rho1*eigen(W1)$values-rho2*eigen(W2)$values))
		lnL <- -(n/2)*log(pi*sigma2)+J-(t(eps)%*%eps)/(2*sigma2)
		return(-lnL)
	}

##########################################################
## MLE for SAR w/ 2 W's [ my.sar2(y,X,W1,W2) ]
## * Maximum likelihood estimation for a SAR model with
## * two spatial weights matrices
## * Takes as an input a vector of outcomes (y),
## * a matrix of predictor variables (X),
## * and the two weights matrices (W1, W2).
## * Output is a list of parameters and summary statistics.
##########################################################


my.sar2 <- function(y,X,W1,W2){
	init <- c(.1,.1,.1,rep(.1,ncol(X)))
	lb <- c(1e-20,-.999,-.999,rep(-1e10,ncol(X)))
	ub <- c(1e20,.999,.999,rep(1e10,ncol(X)))
#pars <- optim(init,y=y,X=X,W=W1,fn=sar.lik,method="L-BFGS-B",lower=lb,upper=ub,hessian=T,control=list(parscale=c(.1,.1,rep(.1,ncol(X)))))
pars <- nlminb(start=init, objective=sar2.lik, gradient = NULL, hessian = sar2.lik, y,X,W1,W2, lower = lb, upper = ub)
hes <- hessian(func=sar2.lik,x=pars$par,y=y,X=X,W1=W1,W2=W2)
s2.hat <- pars$par[1]
rho1.hat <- pars$par[2]
rho2.hat <- pars$par[3]
beta.hat <- pars$par[4:length(pars$par)]
LL <- -pars$objective
aic <- 2*ncol(X)-2*LL
bic <- -2*LL + ncol(X)*nrow(X)
resid <- y-(rho1.hat*W1%*%y+rho2.hat*W2%*%y+X%*%beta.hat)
SSR <- sum(resid^2)
OI <- solve(hes)
se <- sqrt(diag(OI))
coef <- cbind(c(rho1.hat,rho2.hat,beta.hat),se[2:length(se)])
colnames(coef) <- c("Coefficient","Standard Error")
rownames(coef) <- c("Rho1","Rho2",paste("X",1:ncol(X),sep=""))
return(list(pars=pars,coef=coef,s2=s2.hat,LL=LL,AIC=aic,BIC=bic,SSR=SSR))
}

